#pragma once

#include <CVector.h>
#include "CGraph.h"

class C3DField
{
 public:
  //double* Grid; 

  complex<double>* Grid;
  float* dGrid;

  int Width;
  CVector Min, Max, Delta;
  double BoxSize;
  bool PlanFFTInitialized, PlanFFTInvInitialized;

  #ifndef STRIP
  fftw_plan PlanFFT, PlanFFTInv;
  fftw_complex* in;
  fftw_complex* out;
  #endif 



  C3DField();
  ~C3DField();

  void Delete();
  void Initialize(int width, CVector, CVector, double); 
  void InsertDirect(CVector& p, double val); 
  void Clear();
  void Clear(double);
  void toFloatGrid();
  complex<double>& Get(int, int, int); 
  double& GetR(int, int, int); 
  double& GetRBoundary(const int& x, const int& y, const int& z);

  double noval;

  double& GetI(int, int, int); 
  double GetRealVal(CVector& p);
  int GetParticlePosition(CVector& p);
  CVector GetParticlePositionInGrid(CVector& p);
  void Set(int, int, int, double);
  void Set(int, int, int, complex<double>);
#ifndef NO_OPENGL
  void Render();
  
  void toTexture3D();
#endif

  double getAverage() {
    double avg = 0;
    for (int i=1;i<Width*Width*Width;i++)
      avg+=norm(Grid[i]);
//      if (norm(Grid[i])>max) max = norm(Grid[i]);

    avg/=(double)(0.01*Width*Width*Width);
    return avg;
  }

  void Correlation2(string Filename, CGraph& Graph);
  void PowerSpectrum(CGraph& result, CGraph& Growth, int idx);

  void Copy(C3DField& s);
  void Gaussian(double B);
  void Mul(C3DField &other);
  void Scale(double s);
  void PhysicalScale(double s, double rho_c);

  void FFT(C3DField& o);
  void FFTInv(C3DField& o);

  void Differentiate(C3DField& Fx, C3DField& Fy, C3DField& Fz, double cutoff=0);

  #ifndef STRIP
  void SetupFFT(fftw_complex** in, fftw_complex** out, double scale);
  void ResolveFFT(fftw_complex* out, C3DField* );
  #endif STRIP

  void GenerateField(double scale, double D);

};
